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Abstract. On microscopic and mesoscopic scales, plastic flow of crystals is 
characterized by large intrinsic fluctuations. Deformation by crystallographic slip 
occurs in a sequence of intermittent bursts ('slip avalanches') with power-law 
size distribution. In the spatial domain, these avalanches produce characteristic 
deformation patterns in the form of slip lines and slip bands which exhibit long- 
range spatial correlations. We propose a generic continuum model which accounts 
for randomness in the local stress-strain relationships as well as for long-range 
internal stresses that arise from the ensuing plastic strain heterogeneities. The model 
parameters are related to the local dynamics and interactions of lattice dislocations. 
The model explains experimental observations on slip avalanches as well as the 
associated slip and surface pattern morphologies. 
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1. Introduction and Background 

According to the classical paradigm of continuum plasticity, plastic deformation of crystalline 
solids resembles a smooth, quasi laminar flow process. This viewpoint has recently been 

challenged by experimental investigations which have demonstrated that, even in the absence 
of macroscopic deformation instabilities, plastic flow on microscopic and mesoscopic scales 
proceeds in a strongly heterogeneous and intermittent manner. Beyond the microscopic scale, 
where spatio-temporal localization of deformation is a trivial consequence of the discreteness 
of lattice defects (dislocations), slip localization phenomena extend over a wide range of 
mesoscopic scales and may involve the collective dynamics of very large numbers of defects. 
This is well known from the observation of slip lines or slip bands on the surface of deformed 
crystals (for an overview, see Neuhauser, 1984), where slip steps on the surface manifest 
the collective motion of large numbers of dislocations. Recently, it has been shown that the 
spatial arrangement of 'slip events' has a fractal character (Weiss and Marsan, 2003). The 
consequences of a fractal distribution of slip on the surface morphology of plastically deformed 
crystals have been investigated by Zaiser et al. (2004), who demonstrated that the surface of 
plastically deformed Cu samples develops self-afhne roughness over several orders of magnitude 
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Figure 1. Surface profiles of a Cu polycrystal deformed in tension to a plastic strain 
e = 9.6%. Top: AFM profile; bottom: SWLI profile. The x direction is parallel to the 
direction of the tensile axis. After Zaiser et al. (2004) 
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Figure 2. Roughness plots (mean height difference vs. distance along the profile) 
for AFM and SWLI profiles obtained at strains of 9.6 and 17.8 %. The corresponding 
profiles for e = 9.6% are shown in Figure 1. 



in scale. 

Spatio-temporal heterogeneity of slip in spite of macroscopically stable deformation is 
illustrated in Figures 1 and 2. Figure 1 shows surface profiles obtained from a Cu polycrystal 
with an initially flat surface after deformation to 9.6% tensile strain. The statistically self-affine 
nature of the profiles can be inferred by comparing profiles taken by Atomic Force Microscopy 
(AFM, profile length of 25 //m), and by Scanning White-Light Interferometry (SWLI, profile 
length 2 mm). Over several orders of magnitude, the self-affine behavior can be quantitatively 
characterized by a single Hurst exponent H as the average height difference {\y{x) — y{x + L)\) 
between two points on a profile increases as a function of their separation L like with 
H 0.8 (see Figure 2). 

Information about the temporal dynamics of plastic flow can be inferred from acoustic 
emission measurements. Measurements on ice single crystals reported by Weiss and Grasso 
(1997) and Miguel et al. (2001) show temporal intermittency as the acoustic signal consists 
of a random sequence of discrete events ('slip avalanches'). These events exhibit a scale-free 
size distribution, as the probability density to observe events with energy release E decreases 
according to p{E) oc E^^'^ (Figure 3). 

The present paper outlines a simple constitutive model that is capable of capturing spatial 
heterogeneity, temporal intermittency and scale-free behavior in plastic flow of crystals. To this 
end, we adopt concepts from the dynamics of random media by assuming that the deformation 
of a given volume element occurs through a random sequence of hardening and softening 
processes, such that the local flow stress is a fluctuating function of the local strain. Physically, 
these fluctuations result from randomness in the arrangement of the microstructural defects 
(dislocations) that produce the plastic deformation. Because of the crystalline nature of the 
medium, plastic deformation takes place by crystallographic slip, which in the present paper 
is assumed to occur on a single slip system only. Fluctuations in the local flow stress lead to 
shear strain fluctuations. These, in turn, give rise to long-range stress redistribution which 
can be expressed in terms of an elastic Green's function. Furthermore, dislocation-dislocation 
correlations give rise to a local back stress which can be expressed in terms of a second-order 
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Figure 3. Distribution of energy releases in acoustic emission during creep 
deformation of ice single crystals; temperature T = 263 K, resolved shear stresses 
on the basal plane as indicated in the inset. After Miguel et al. (2001). 



gradient of strain. These ideas are detailed in Section 2 where a stochastic plasticity model 

is presented. In Section 3, the numerical procedure employed to implement the model is 
briefly outlined. In Section 4, we discuss results obtained from the numerical simulations, and 
demonstrate that the model correctly accounts for the observed statistics of 'slip avalanches', 
the formation and arrangement of slip lines, and the corresponding deformation-induced 
surface morphology. 

2. A Crystal Plasticity Model with Disorder 

2.1. Basic structure of the model 

We consider the simplest possible model for crystal plasticity with plastic deformation 
occurring by crystallographic slip on a single active slip system. Hence, the plastic distortion 
tensor is given by (3^{r) = j{r)n s where 7 is the shear strain on the slip system, and n 
and s are unit vectors pointing in the direction of the slip plane normal and the slip direction, 
respectively. In the following we assume without loss of generality that the slip direction 
corresponds to the x direction of a Cartesian coordinate system and the slip plane is the xz 
plane. The driving force for plastic flow is the shear stress (Txy{r) =: r(r) acting in this slip 
system. We consider rate-independent elastic-plastic behavior. Hence, the stress acting in a 
volume element located at r must fulfill the inequality 

r(r) < Tf(r,7(r)) . (1) 

Here Tf is the local flow stress, which depends on internal state of the volume element under 
consideration. It is in general a function of the local strain 7(r) and may also depend explicitly 
on the space coordinate r. The stress r(r) = Text+Tint('") acting from outside on the considered 
volume element is a sum of external and long-range internal stresses. The external stress, which 
acts as an external driving force on the system, is assumed space independent over the region 
of interest. The space-dependent internal stress field Ti^tif) is a functional of the (in general 
spatially non-homogeneous) strain field 7(r); it vanishes for a homogeneous deformation state. 
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If the inequality (1) is violated, the local strain 7(r) increases quasi-instantaneously until 
Eq. (1) is again satisfied. Impossibility to satisfy this equation implies ductile failure of the 
system. 

2.2. Evaluation of long-range internal stresses 

We calculate the internal stresses in an infinite three-dimensional body with an arbitrary 
plastic distortion field 13^ [r). The external stress is assumed to be zero. (A non-zero external 
stress simply adds to the internal stresses.) We start out from the elastic equilibrium equation 
for the components aij of the stress tensor: 

-djaij = fi , (2) 

where fi are the body forces and sums are performed over repeated indices. The above equation 
can be rewritten in terms of the components Ui of the elastic displacement vector as 

—djCiikidkUi = fi, (3) 

where Cijki are components of Hooke's tensor. The solution of Equation (3) has the form 

Uiir) = J r,fc(r-r')/fe(r')dV (4) 

where the Fourier transform of the elastic Green's tensor Tik(r) is rjfc(k) = [—Cijkikjki]^^ . 

We now first consider the particular problem of a plastically deformed 'inclusion' where 
the plastic distortion has a constant value (3^ over a certain volume V and is zero elsewhere. 
This inclusion problem is solved as follows: The volume V is first cut out of the surrounding 
matrix and deformed plastically in order to produce a stress- free strain /?p. To re-insert it into 
the matrix, interface tractions arc applied, such that the original shape is restored. According 
to Equation (2) these interface tractions are 

fi = CijkiPlidjHvir), (5) 

where the function Hy{r) is equal to unity within V and zero elsewhere. The volume is then 
placed in its original position and relaxed. Elastic relaxation proceeds until the tractions 
produced by the relaxation strain balance those given by Equation (5). The Fourier 
transform of the corresponding displacement field is 

<(k) = ^^ik{k)CklmnP^nklHv{^<.) (6) 

The total elastic distortion is then the sum of the relaxation strain and the initial distortion 
— /3?--ffy(r) applied to "restore" the original shape before relaxation. The associated total 
stress reads crij(r) = CijkiiPli ~ (^kl^v{T^)) and its Fourier transform is 



The angular average 



(k) = —Cijim (^kmkoi'ln{k)Cnopq + g^lp^mq^ f^pq^V ■ (7) 



~Q 1 f ~ 1 

^ijlm ~ ^ / C^ijno[kokq^npO^)Cpqlm + Q^ln^mold^ (8) 

does not depend on the modulus k of the wavevector since Fnp scales like k~'^. The stress 
hence can be written as 

ai,{k) = [-fo,„ + t*{k/k)]Pf^Hv{k) , (9) 
where the second term is simply defined by f * (k/fc) = Cijno[kokqfnp{k)Cpqim+lSin6mo\-i'^jim- 
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The above procedure can be straightforwardly generahzed for an arbitrary distribution 
of the plastic distortion (r) by considering each volume element as a separate "inclusion" . 
Equation (9) becomes 

ai,(k) = [-f + f *(k/fc)]/3P„(k) (10) 

and, in real space, 

,(r) = -f°,„/3L(r) + / r*(r - r')C(r')dV' . (11) 

The non-local kernel r*(r — r') is the inverse Fourier transform of r*(k/A;); it scales like 
and it can be shown that it has zero angular average. What has been done so far holds for 
strain fields that go to zero at infinite distances. If the asymptotic value of the plastic strain 
assumes a non-zero value we have to add the corresponding stress-free strain as follows: 

(k) = - + I r*(r - r')/3L(r')rf'r'. (12) 

For a plastic distortion field which has the average value (/S^), the asymptotic value /3P'°° is 
replaced by the average (/?p) since the fluctuation contributions average out if integrated over 
the infinite contour. Hence, the internal stress can be envisaged as the sum of a mean-field 
contribution and a non-local term with a kernel of zero angular average. For the purpose of 
a depinning theory built with this type of non-local elastic interaction, one may note that in 
Fourier space the kernel scales like with a = 0, i.e., mean-field theory is expected to be 
valid in all dimensions. 

For later use we give explicit expressions for the case where the plastic strain is determined 
by slip on a single slip system, = j{r)ey (g) e^, and the shear strain 7 depends on the x and 
y coordinates only. (Such a quasi-two-dimensional model corresponds to a system of straight 
parallel edge dislocations, cf. below.) In this case, the internal shear stress Tint = CTxy is in 
Fourier space given by 

and in real space by 



(^int(r-)) = / lir') 



8{x - x'fiy - 



27r(l -z/) J " ' [{r-r'f (r 



d r 



+ 4(l^[(7)-7(r)]. (14) 

Here, G is the shear modulus and v the Poisson's ratio. Prom these expressions, two points may 
be noted: (i) The elastic kernel is not positively definite in real space, (ii) There exist certain 

space-dependent strain fluctuations (fluctuations with wavevectors in the x and y directions) 
which do not give rise to any long-range internal stresses. The implications of this will be 
discussed below. 



2.3. Dislocation-related stresses and internal- stress fluctuations 

We now discuss the physical origin of the 'flow stress' Xf (r, 7). On the microscopic scale, plastic 
flow of a crystal is brought about by the motion of linear lattice defects (dislocations). The 

flow stress of a small mcsoscopic volume (a volume containing multiple dislocations/dislocation 
segments) corresponds to the stress required to move dislocations through the stress 'landscape' 
within this volume. In the absence of other defects, this 'landscape' is created by the 
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dislocations themselves. The resulting flow stress T{ = St + Tp can be envisaged as a sum 
of two contributions [for a more detailed discussion, see Zaiser and Seeger (2002), Groma, 
Csikor and Zaiser (2003), Zaiser and Aifantis (2005)]: 

(i) A spatially fluctuating stress 5r(r) which depends on the positions of the individual 
dislocations. The average of the fluctuating stress is zero, and the characteristic 'amplitude' 
of the fluctuations is given by 

{6T{rf) = K^G^h^pir) , (15) 

where K'^ = ln(^/6)/[87r(l — u)"^] is a numerical constant of the order of unity which depends 
on the elastic properties of the crystal lattice, and logarithmically on the characteristic range 
^ « of dislocation-dislocation correlations. The spatial correlation function of the 

fluctuating internal stress field is given by 

{5T{r)5T{r + r')) = {5T{rf)h{r'/0 , (16) 

where h{r' /E) is a short-range correlation function with characteristic range ^ « (Zaiser 
and Seeger, 2002). 

Dislocation glide increases the local strain and modifies the fluctuating stresses within 
the surrounding mesoscopic volume element. In the absence of detailed information about the 
individual dislocation positions, we take this evolution of the stress 'landscape' into account 
by envisaging the fluctuating stress as a random function not only of space but also of strain, 
with the correlation function 

(<5r(r, 7)^V(r + r', 7 + 7')) = {5r{rf)h{r' /Ooil llcorr) , (17) 

where g is another short-range function. The 'correlation strain' 7corr can be estimated as 
the local strain produced when all dislocations in a volume element move by the distance ^, 

Tcorr = ~ b^. 

(ii) If small groups of dislocations move in a correlated manner, their mutual interactions 
tend to homogenize deformation in the slip plane. This can be described by a 'pile-up stress' 
Tp{r) which can be approximated by a second-order gradient of the strain according to 

DG DG 
Tp(r) = —Ixx id = 2) , Tp(r) = —[yxx + Izz] (c? = 3) . (18) 

Here, p is the total dislocation density and D a constant of the order of unity [see Groma, 
Csikor and Zaiser (2003)]. In Eq. (18), d € [2,3] is the dimensionality of the model. The first 
expression (d = 2) refers to the special case where deformation is due to the motion of straight 
edge dislocations, such that the shear strain 7 depends on the x and y coordinates only. 



2.4- Strain hardening 

During deformation, the dislocation density often increases with increasing strain, leading to 
an increase of the stress required to sustain plastic flow (strain hardening). In our model we 

may take this into account in two manners: 

(a) We explicitly account for the strain evolution of the dislocation density p. In the 
simplest case, this is given by [see e.g. Zaiser (1998)] 

where /3 0.02 is a small non-dimensional parameter. It follows that the dislocation density is 
approximately constant over the 'correlation strain' 7corr ~ but increases slowly over larger 
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strain intervals. This leads to a corresponding increase of the amplitude of the fluctuating stress 
5T(r,7), which scales in proportion with the local strain 7(r). 

(b) In phenomenological hardening theories, hardening is often described by a strain 
dependent 'back stress' which is subtracted from the locally acting stress. In the simplest 
case of linear strain hardening, this stress is given by vit, = —87 where B is the hardening 
coefficient. If we choose this description of hardening, we keep the amplitude of the fluctuating 
stresses constant but subtract a slowly increasing back stress which is proportional to the local 
strain. 

Actually, both descriptions yield very similar results. As discussed in the next paragraph, 
our equations have the structure of an clastic manifold depinning model, with the fluctuating 
stress St acting as a random pinning field. Since the pinning stress in such models is governed 
by the negative maxima of the pinning field, a linear increase of the amplitude of the pinning 
field with increasing strain, or the subtraction of a strain-dependent back stress, both have 
the same consequence: The stress required to overcome these maxima increases linearly as a 
function of the local strain. 



2. 5. Plastic flow and elastic manifold depinning 

Collecting all stress contributions discussed in Sections above, we find that plastic flow occurs 
(7 increases quasi-instantancously) as soon as the inequality 

DG 

Text + Tint(r) + — [^ixx + lyy] + (^T(r, 7) < (20) 

is violated. The asymptotic behavior of this model for positive external stresses corresponds 
to that of the equation of motion 

1 DG 

— dtjir) = Text + Tint{r) H [jxx + lyy] + ST{r, 7) . (21) 

B p 

in the rate-independent limit B oo. Formally, Eq. (21) can be envisaged as describing 
the overdamped dynamics of an elastic manifold with coordinates (7, r) which moves in 
the 7-direction through a random medium. In the two-dimensional case where 7 = 7(x,y) 
depends only on the x and y coordinates (d = 2), our model may be considered a continuum 
approximation of a quasi-two-dimensional system of straight parallel dislocations. If the strain 
depends on all three spatial coordinates, and the model mimics the large-scale behavior of a 
system of three-dimensionally curved dislocations. In either case, the model can be formally 
envisaged as describing the motion of a d-dimensional elastic manifold through a disordered 
d + 1-dimensional medium which exerts a fluctuating pinning force St. Due to the inflnite 
range of the interaction kernel governing the internal stress rint(*'), one expects the second- 
order gradient term in Eqs. (20,21) to be irrelevant for the large-scale behavior, and the model 
exhibits mean-field behavior irrespective of the dimensionality of the 'manifold'. 

Two differences may be noted with respect to conventional models of elastic manifold 
depinning: (i) If we account for hardening, cither the mean value or the amplitude of the 
pinning field increases as the manifold advances, i.e., we are dealing with a net 'uphill' motion 
of the manifold. The consequences of this will be discussed in Section 4.2. (ii) As can be 
directly seen from Eq. (14), the elastic kernel mediating the long-range interaction is not 
positively definite. Hence, the 'no passing'-theorem (Middleton 1992) need not hold and the 
existence of a unique depinning threshold is not a priori guaranteed. In spite of these caveats, 
numerical investigation of the model shows depinning-like behavior, as demonstrated in Section 
4. 
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3. Numerical Implementation 

For numerical implementation of the model we focus on the two-dimensional case. (Because 
of the mean-field nature of the long-range interactions, a generalization to three dimension is 
not expected to make any significant difference). We use a lattice automaton model where we 
discretize space in terms of a two-dimensional array of cells, each cell representing a mesoscopic 
volume element whose edge length / « ^ is of the order of the correlation length of the stress 
fluctuations. Accordingly, the fluctuating stresses in different volume elements are assumed as 
statistically independent random variables. In a similar spirit, the local strains are discretized 
in units of the correlation strain 7corr ~ i>\/Pi and the fluctuations at difl"crcnt strain steps 
are considered statistically independent. A bulk system is mimicked by imposing periodic 
boundary conditions in both directions. The total dislocation density p and, hence, the 
statistical properties of the fluctuating internal stress field are assumed spatially homogeneous. 
A non-dimensional and dislocation density independent formulation is achieved by scaling all 
lengths in proportion with ^ « all stresses in proportion with KGb^yp, and all strains 

in proportion with 6/^ « i>^/p. 

The procedure for implementing the random function (5t(7, r) in Eq. (20) is here discussed 
in the absence of hardening. To each cell we assign initially a negative random value of 5t 
drawn from the negative half of a Gaussian distribution with unit variance [Eq. (3)]. This 
ensures that there is no plastic flow at zero external stress; in physical terms, it implies that 
initially all dislocations are trapped in conflgurations where they are held up by negative (back) 
stresses. Whenever the strain within a cell increases, we assign a new (positive or negative) 
random value of 5t to this cell, which we draw again from a Gaussian distribution of unit 
variance. We note in passing that the Gaussian shape of the probability distribution of 6t is 
not crucial; other distributions with finite first and second moments produce similar results. 

In each simulation we increase the applied stress Text from zero in small increments Arext • 
After a stress increment we check for all volume elements whether Eq. (20) is violated. 
In all volume elements where this happens, we increase the local strains by a unit amount 
and assign new values of 6t. After this simultaneous update procedure is completed , we 
compute the internal stress field (rint(T')) corresponding to the new strain pattern by using 
the periodically continued elastic Green's function of a two-dimensional isotropic medium, Eq. 
(14), and evaluate the gradient-dependent stress contribution using the discrete second-order 
gradient. Once the new internal stresses are computed, we check again for all volume elements 
whether Eq. (20) is now satisfied and, if not, increase the local strains etc.. The process 
is repeated until the local stress everywhere falls below the local yield stress, or until the 
average strain within the system exceeds a prescribed maximum value at which we terminate 
our simulation. If the system settles into an equilibrium configuration before reaching the 
terminal strain, we increase the external stress further by Arext) and so on. This procedure 
implements an adiabatically slow drive, provided that the external stress increment Arext is 
chosen small enough (the practical criterion being simply that the results do not change if the 
stress increment is further decreased). 

4. Results and Discussion 

4-1- Stress-strain characteristics and avalanche dynamics in the absence of hardening 

We first study the model in the absence of hardening. Then, the 'pinning stress' St has 
stationary stochastic properties. A stress-strain graph obtained for this case is illustrated in 
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Figure 4. The calculation has been done for the following parameters typical for the onset 
of deformation of Cu single crystals: G = 40000 MPa, u = 0.3, 6 = 2.5 x 10~i°m, p = 10^^ 
m^^, and ^ = = 10~^m. With these parameters, ~ 1, while the parameter D in 

Eq. (18) was chosen as D = 0.1. We note that, in the context of the present model, this 
latter term mainly serves to break the symmetry which exists in the elastic kernel given by 
Eq. (14) between the directions in and out of the slip plane. The numerical value of D is of 
little influence on our results unless very high D values are assumed. 

From Figure 4, the following observations can be made: (i) In spite of the fact that 
the averages of St over 7 at fixed r (or over r at fixed 7) are zero, a finite stress (yield 
stress) of Tc ~ 0.4 (in scaled coordinates) is required to sustain deformation. This is because 
at low applied stresses the system becomes pinned in configurations where the fluctuating 
stress 6t is negative in most volume elements, thereby creating a non-zero average back stress, 
(ii) The yield stress is approached asymptotically through a 'microplastic' region where such 
metastable configurations are gradually eliminated, (iii) The increase of plastic strain with 
increasing stress occurs in discrete 'slip avalanches' of varying size. These avalanches arc visible 
as steps on the stress-strain curves, which assume a staircase-like shape. The intervals between 
the larger avalanches divide into avalanches of smaller size, and the characteristic avalanche 
size diverges as one approaches the yield stress. In dimensional units, the scaled stress and 
strain ranges visible in Figure 4 correspond to < Text < 5 MPa and < (7) < 0.5%, 
respectively. The figure can, hence, be thought of as describing the microplastic region and 
yielding transition at the onset of easy-glide deformation of a fee single crystal. We also note 
that the dimensional yield stress obeys Taylor's well-established scaling relation Tc = aGby/p. 
The numerical value of a ~ 0.4 obtained from the simulation is within the characteristic range 
0.2 < a < 0.5 observed in experiment. 

We first study the average behavior of our model. By averaging the stress-strain graphs 
over a large number of simulations, we obtain a (nearly) smooth stress-strain relationship 
(Figure 5). The average strain diverges as the stress approaches the yield stress Tc- A 
semi-logarithmic plot of 7 vs. Tc — Text reveals that this divergence is logarithmic in nature. 



0.5 




Figure 4. Stress-strain curve as obtained from simulation of a system with 128 x 128 
sites; dashed line: critical stress t^, Insert: detail of the same stress-strain graph. 
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Figure 5. Ensemble-averaged stress-strain graph as obtained by averaging over 60 
simulations of systems with 128 x 128 sites. Strain is plotted against distance from 
the critical stress. 
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Figure 6. Probability distributions of slip avalanche sizes (probability density p{AE) 
vs. energy release AE) as obtained from an ensemble of systems of size 256x256. Left: 

distributions corresponding to different stresses: Right: vmiversal distribution obtained 
by re-scaling AE AE{1 — Text/Tc)^; full hue: theoretical curve (Eq. 23). 



Accordingly, the stress susceptibility x = dj/dText of the plastic strain diverges according to 

X OC (Tc - Text)"^ ,0^1. (22) 

We now investigate in more detail the 'strain bursts' that manifest themselves as stepwise 
strain increments in the individual stress-strain graphs. To this end, we record the strain 
increments (if any) which occur after each stress increment Axext- Since the total strain 
increment due to a slip avalanche depends on the system size Vs (a given local slip avalanche 
produces a smaller total strain increment A7 in a larger system) , we characterize the avalanche 
sizes in terms of the dissipated energy AE = rextA7Vs which is proportional to the strain 
increment but does not depend on system size. We then determine the probability distribution 
of avalanche sizes within different external stress intervals, corresponding to different distances 
from the critical stress Tc- Results are given in Figure 6. 
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The distributions exhibit a power-law decay p{AE) oc AE~'^ which is truncated at a 
characteristic avalanche size AEc- As the stress approaches the critical stress Tc, this upper 
limit and hence the average avalanche size diverge. Distributions obtained at different applied 
stresses Text can be collapsed into a single universal distribution if the energy releases are 
re-scaled by a factor (1 — Text/Tc)^^'^ where a fa 0.5 (Figure 6, right). The results for 
the susceptibility and the avalanche statistics are not independent since the susceptibility 
is proportional to the average avalanche size. We therefore expect it to diverge like x 
(tc - t)~^'^~'^^/'^. With K 1.4 and a ^ 0.5 we find = (2 — K)/cr 1 in agreement with Eq. 
(22). The exponents k, 9 and a can be understood in terms of mean-field depinning [see e.g. 
Fisher (1998) or Zapperi et al. (1998)], as the theoretical values 9 = l,a = 0.5 and k = 1.5 are 
in good agreement with the present simulation results. The theoretically predicted avalanche 
size distribution for mean-field depinning is 

2' 



p{AE) oc A^;-^-^ exp 



/ AE 

[ae. 



(23) 



which again compares well with the results of the simulations (Figure 6) . 



4-2. Influence of hardening 

To introduce hardening with non-dimensional hardening coefficient into our model, we use 
two approaches: (i) We increase the amplitude of the internal-stress fluctuations (i.e., the 
width of the probability distribution of St) in proportion with the local strain by an amount 
of Q'y/Tc{Q = 0); (ii) We leave the width of the distribution unchanged but shift the mean 
value by an amount 07 that is again proportional to the local strain. Both methods result in 
stress-strain graphs similar to those depicted in Figure 7: The asymptotic behavior changes, as 
the horizontal tangent of the stress-strain graph is replaced by a tangent with slope 0. To give 
again a 'feeling' for the orders of magnitude involved, we note that the simulated hardening 
rates correspond, with the parameters for Cu given above, to dimensional hardening rates of 




4 8 12 16 20 24 28 32 36 40 

plastic Strain [bp^'^] 

Figure 7. Stress-strain curves as in Figure 4 for hardening coefBcients O = 0.001 and 
6 = 0.004; dashed lines: curves t = Tc + 67 with = 0.44. 
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Figure 8. Probability distributions of slip avalanche sizes (probability density p{AE) 

vs. energy release AE') as obtained from an ensemble of systems of size 256x256. Left: 
distributions corresponding to different hardening rates; Right: universal distribution 
obtained by re-scaling AE — > QAE 



40 and 200MPa which are in the range of what is typically observed in hardening stages I 
and II of Cu single crystals. 

It is seen from Figure 7 that the influence of hardening does not change the stepwise 
morphology of the stress-strain curves but eliminates the largest strain bursts which, in a non- 
hardening system, occur at stresses close to the critical stress. To quantify this idea, we study 
the avalanche statistics (evaluated for stresses above the critical stress of the non-hardening 
system) for different hardening rates. Figure 8 shows that the hardening does not change the 
power-law regime of the burst size distribution but decreases the cut-off value, which turns 
out to be inversely proportional to the hardening rate. 

This can be easily understood from a scaling argument. The burst energy is proportional 
to the strain produced in a burst, and with an avalanche exponent oi k ^ 1.5 it follows that the 

1/2 

average strain A7 follows the relation (A7) oc {A.E) oc AE'max. Due to hardening, each burst 
raises the critical stress by a small amount Atc oc Q{Aj), which implies that the external 
stress (which cannot instantaneously follow) lags behind the current critical stress by a similar 
amount. Finally, the maximum burst size follows the scaling relationship AE^^x oc Ar~^. 
Combining these relations, it follows directly that AE'max oc 1/6. 

4-3. Slip pattern and surface morphology 

As shown in the previous sections, mean-field depinning provides a framework for 
understanding the size distribution of slip avalanches and the behavior of the stress-strain 
curves as one crosses the microplastic region and approaches the yield stress, as well as for 
assessing the influence of strain hardening. However, mean-field considerations can in principle 
not account for the observed anisotropic spatial distribution of slip, the formation of slip lines, 
and the development of a self-affinc surface morphology. However, our simulations indicate 
that, owing to the particular properties of the elastic kernel, such features naturally emerge 
from the dynamics of our model. 

Figure 9 shows a grey-scale representation of the slip pattern ^{x, y) that emerges in a 
typical simulation. One observes a 'slip-line'-pattern with pronounced correlations along the 
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Figure 9. Strain pattern obtained after simulation of a system of size 256 x 256 to 
an average strain of 20by/p (slip direction from left to right); parameters as in Figure 
4; greyscale: local strain in units of 



direction of slip (x direction in Figure 9). Such stratified slip patterns are commonly observed 
in deformation of single-slip oriented single crystals. They are considered to be a natural 
consequence of the directed glide of lattice dislocations along the slip planes. The present 
model produces 'slip-line' patterns without explicitly accounting for the motion of lattice 
dislocations. This is can be understood by noting that according to Eq. (13)) heterogeneities 
of slip in the direction of the slip plane normal (the y direction) do not give rise to any long- 
range stresses. Heterogeneities in the slip direction (the x direction), on the other hand, are 
penalized by the gradient-dependent stress contribution in Eq. (20) which tends to homogenize 
deformation in the slip direction. Because of these properties of the interactions, the 'manifold' 
7(x, y) is prone to fall apart in the y direction, leading to the striated deformation patterns 
seen in Figure 9. 

We now evaluate the surface profiles which correspond to the simulated slip patterns. We 
consider a surface running along the plane x = 0. The displacement u in the direction normal 
to the surface fulfils the relationship 

^ = j{x = 0,y). (24) 

This can be understood as the sum of a displacement due to the average strain (7) , which leads 

to a rigid rotation of the surface, plus a random displacement due to the strain fluctuation 
5j{x = 0, y) = 7(x = 0,y) — (7). The surface profile may be obtained by direct integration of 
the strain fluctuation, 

%)= r[^{x = 0,y')-{^)]dy' . (25) 
Jo 

A surface profile obtained in this manner is shown in Figure 10. 

To assess simulated profiles in terms of their roughness exponent, we plot the mean square 
height difference {\u{x) — u{x + between two points along the profile as a function of their 
separation I. For a self-affine profile we expect {\u{x) — u{x + l)\) oc l^ where H is the Hurst 
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Figure 10. Surface profile obtained from a system of size 512 x 512 after simulation 

to an average strain of 7 = 20by/p. 



0) 

o 
a5 



w 5 

^ 10 



2 

CO 

w 
c 

CO 
CD 



10^ 



10" 



H=0.7 




Average strain 6bp 

Average strain 2Qbp' 



10 100 
Horizontal spacing [^] 



1000 



Figure 11. Mean square height difference vs. horizontal distance for surface profiles 
obtained from a system of size 512 x 512 after simulation to average strains of Qb^/p 
and 206y^; each graph has been averaged over 10 simulated profiles. 



exponent. Figure 11 demonstrates that our simulated profiles can be characterized by a Hurst 
exponent H 0.7. Increasing the total strain leads to an increase in the absolute magnitude of 
surface height variations but docs not change the value of H. Both the numerical value of the 
Hurst exponent and its strain independence compare well with the experimental observations 
quoted in the introduction. 



5. Conclusions 



Our simulations demonstrate that it is possible to capture essential aspects of slip avalanches in 
terms of a simple continuum model. By combining a fluctuating local stress-strain relationship 
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with a strain-gradient dependent stress contribution and long-range stresses mediated by the 
elastic Green's function, the model accounts both for the observed spatial heterogeneities of 
plastic deformation (slip lines, self-afiine surface roughness) and the emergence of 'bursts' of 
plastic activity with a power-law size distribution. 

The present model has a fairly generic formal structure, and many of our qualitative 
results do not depend critically on the particular physical assumptions made. For instance, it 
is not crucial to relate the fluctuating flow stress to the stress fluctuations associated with a 
homogeneous distribution of dislocations as done in Section 2.3. Concepts of randomness 
and disorder also apply to flow stress fluctuations arising from other types of random 
microstructurc such as irregular patterns of dislocation cells, or from the statistical activation 
and exhaustion/blocking of dislocation sources. However, there are still several crucial 
questions that need to be studied: (i) Is it possible to formulate a dislocation-based model 
which exhibits the same behavior as the present phenomenological continuum model? (ii) How 
does the behavior of the model change if we allow for slip on multiple slip systems, such that 
the tensorial nature of the plastic strain becomes important? (iii) What is the influence of 
surface boundary conditions on the slip patterns - are the simulated 'surface profiles' truly 
equivalent to those emerging on a free surface, where the stresses are by necessity different 
from the bulk? 

In general terms, or results indicate that fluctuation phenomena in plastic flow of 
crystalline solids can be envisaged as critical phenomena in driven non-equilibrium systems. 
The theory of elastic interface depinning provides a theoretical framework which covers a 
wide range of phenomena including avalanche dynamics in the motion of magnetic domain 
walls (Barkhausen noise, Zapperi et al. 1998), the dynamics of earthquakes or the motion of 
cracks (Fisher, 1998). The present model adopts this framework to problems of plastic flow. 
An analogous effort has been made in relation with the plasticity of amorphous materials by 
Baret et al. (2002). Peculiarities arise in the present case from the particular deformation 
geometry: In a crystalline solid deforming in single slip, where deformation is restricted to 
simple shear occurring along a single set of planes, the elastic kernel which mediates long-range 
elastic interactions is not positively definite. The consequences of this observation as well as 
the reasons for the emergence of a power-law decay of correlations in the strain pattern (which 
has, to the knowledge of the authors, not been observed in other depinning-type models) 
remain to be assessed in future work. 
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